#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 19 11:35:50 2023

@author: jcfq2
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import path



def inpolygon(data_x, data_y, mask_x, mask_y):
    shape = data_x.shape
    data_x = data_x.reshape(-1)
    # print(data_x.shape)
    data_y = data_y.reshape(-1)
    mask_x = mask_x.reshape(-1)
    mask_y = mask_y.reshape(-1)
    q = [(data_x[i], data_y[i]) for i in range(data_x.shape[0])]
    p = path.Path([(mask_x[i], mask_y[i]) for i in range(mask_x.shape[0])])
    return p.contains_points(q).reshape(shape)



# def maskeqzero(data_x, data_y, data_v, mask_x, mask_y):
    
#     inpoly=(inpolygon(data_x, data_y, mask_x, mask_y))
#     data_v[inpoly]=0.
#     return data_v

def maskeqvalue(data_x, data_y, data_v, mask_x, mask_y,setvalue=0.0,invert=False):

    inpoly=inpolygon(data_x, data_y, mask_x, mask_y)
    data_v2=np.zeros_like(data_v) #maybe piped out mask_v with zeros? so added this
    data_v2[:]=data_v[:]

    if invert == True: inpoly=np.invert(inpoly)
    data_v[inpoly]=setvalue
    return data_v


def masksubarray(data_x, data_y, data_v, mask_x, mask_y,invert=False):
    
    inpoly=inpolygon(data_x, data_y, mask_x, mask_y)
    
    if invert == True: inpoly=np.invert(inpoly)
    data_x2=data_x[inpoly]
    data_y2=data_y[inpoly]
    data_v2=data_v[inpoly]
    return data_v2,data_x2,data_y2




# polygon surrounding points - draws a polygon around some stuff

import pylab as p

def _angle_to_point(point, centre):
    '''calculate angle in 2-D between points and x axis'''
    delta = point - centre
    res = np.arctan(delta[1] / delta[0])
    if delta[0] < 0:
        res += np.pi
    return res

def _draw_triangle(p1, p2, p3, **kwargs):
    tmp = np.vstack((p1,p2,p3))
    x,y = [x[0] for x in zip(tmp.transpose())]
    p.fill(x,y, **kwargs)

def area_of_triangle(p1, p2, p3):
    '''calculate area of any triangle given co-ordinates of the corners'''
    return np.linalg.norm(np.cross((p2 - p1), (p3 - p1)))/2.


def convex_hull(points, smidgen=0.0075):
    '''
    Calculate subset of points that make a convex hull around points
    Recursively eliminates points that lie inside two neighbouring points until only convex hull is remaining.

    :Parameters:
    points : ndarray (2 x m)
    array of points for which to find hull
    graphic : bool
    use pylab to show progress?
    smidgen : float
    offset for graphic number labels - useful values depend on your data range

    :Returns:
    hull_points : ndarray (2 x n)
    convex hull surrounding points
    '''

    n_pts = points.shape[1]
    assert(n_pts > 5)
    centre = points.mean(1)
    angles = np.apply_along_axis(_angle_to_point, 0, points, centre)
    pts_ord = points[:,angles.argsort()]
    pts = [x[0] for x in zip(pts_ord.transpose())]
    prev_pts = len(pts) + 1
    k = 0
    while prev_pts > n_pts:
        prev_pts = n_pts
        n_pts = len(pts)
        i = -2
        while i < (n_pts - 2):
            Aij = area_of_triangle(centre, pts[i],     pts[(i + 1) % n_pts])
            Ajk = area_of_triangle(centre, pts[(i + 1) % n_pts], \
                                   pts[(i + 2) % n_pts])
            Aik = area_of_triangle(centre, pts[i],     pts[(i + 2) % n_pts])
            if Aij + Ajk < Aik:
                del pts[i+1]
            i += 1
            n_pts = len(pts)
        k += 1
    return np.asarray(pts)


def polygon_surrounding_points(x,y,padding=1.0,xsize=100,smoothing=0):
    # padding scales from zero, so if you are offset, it breaks

    import scipy.interpolate as interpolate

#    fig = p.figure(figsize=(10,10))


    points = np.ndarray((2,len(x)))
 
    
    points[0,:],points[1,:] = x,y



    padding = padding
    hull_pts = padding*convex_hull(points)

    print(hull_pts.shape)

    p.plot(x,y,'ko')

    x,y = [],[]
    convex = padding*hull_pts
    for point in convex:
        x.append(point[0])
        y.append(point[1])
    x.append(convex[0][0])
    y.append(convex[0][1])

    x,y = np.array(x),np.array(y)

    
    p.scatter(x, y,c='r',linewidth=2)

    
    xsize=xsize
    tck, u = interpolate.splprep([x, y], s=smoothing, per=True)
    x2, y2 = interpolate.splev(np.linspace(0, 1, xsize), tck)

    return x2,y2






# theta = 2*np.pi*np.random.rand(1000)
# r = np.random.rand(1000)**0.5
# x,y = r*p.cos(theta),r*p.sin(theta)

# x2,y2 = polygon_surrounding_points(x,y,padding=1.01,xsize=500)

# p.plot(x2, y2,'r--',linewidth=2)

# p.show()




'''
xpos = np.array([[0.1,0.5,0.9,0.2,0.4,0.5,0.5,0.9,0.6,0.8,0.7,0.2],[0.15,0.55,0.95,0.25,0.45,0.55,0.55,0.95,0.65,0.85,0.75,0.25]])
ypos = np.array([[0.4,0.6,0.9,0.7,0.3,0.8,0.2,0.4,0.4,0.6,0.2,0.6],[0.4,0.6,0.9,0.7,0.3,0.8,0.2,0.4,0.4,0.6,0.2,0.6]])
colorar=np.zeros_like(xpos)+1
for iy, ix in np.ndindex(colorar.shape):
    colorar[iy,ix] = (np.random.rand()*100)+100


mask_x = np.array([0.05,0.4,0.7,0.9,0.05])
mask_y = np.array([0.4,0.7,0.8,0.2,0.4])


colorar=maskeqzero(xpos, ypos, colorar, mask_x, mask_y)


plt.plot(mask_x,mask_y)
plt.scatter(xpos,ypos)

plt.scatter(xpos,ypos,c=colorar)
'''


